/*	UNIVERSIDADE DE SÃO PAULO
	INSTITUTO DE CIÊNCIAS MATEMÁTICAS E DE COMPUTAÇÃO
	
	TRABALHO DE RECUPERAÇÃO - SSC0143 PROGRAMAÇÃO CONCORRENTE
	DESAFIO ERAD_SP 2014 - LANÇAMENTO BALÍSTICO
	profº: JÚLIO CEZAR ESTRELLA

	ALUNO: VITOR BARBOSA NOVELLO
	NºUSP: 7656425

	MAIS INFORMAÇÕES EM
	http://eventos.ufabc.edu.br/erad/problema.pdf

	SÃO CARLOS
	FEVEREIRO/2015
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <pthread.h>

#define pi 3.14159265358979323846
#define NUM_ANGLES 159
#define NUM_THREADS 4

// Estrutura da parâmetro das threads
struct parameter{
	double *s_angle;
	int thread_id;
};

const double g0 = 9.80665;	// Gravidade inicial (nível do mar)
const double R = 6335437;	// Raio da Terra
const double p0 = 1.2922;	// Densidade inicial do ar (nível do mar)
const double e = 0.01;		// Erro máximo tolerado
const double h = 0.0001;	// Passo de integracao para o método Runge-Kutta

// Entradas
double m0;	// Massa total do foguete
double mp;	// Massa do combustível
double tp;	// Tempo de queima do combustível
double Fp;	// Força do empuxo do combustível
double a;	// Ângulo entre vetor empuxo e a tangente local
double Cd;	// Coeficiente de arrasto aerodinâmico	
double r;	// Raio do círculo da secção transversal do foguete
double Sf;	// Alcance horizontal final

// Funções ajustes de ...
double gravity(double);					// ...gravidade
double density(double);					// ...densidade do ar
double mass(double);					// ...massa do foguete
double force(double);					// ...força de empuxo
double drag(double, double, double);	// ...força de arrasto

double ODE_f1(double, double, double, double);	// Função da EDO 1 (aceleração)
double ODE_f2(double, double, double, double, double);	// Função da EDO 2 (rotação)
double ODE_f3(double, double);	// Função da EDO 3 (alcance horizontal)
double ODE_f4(double, double);	// Função da EDO 4 (alcance vertical)

// Função das threads
void* calculate(void*);

// Método Runge-Kutta para calculo de alcance, dado um ângulo
double runge_kutta(double);


int main(int argc, char* argv[]){
	
	// Variáveis auxiliares
	char mode;
	int i, found=0;

	double s_angle[NUM_ANGLES];		// Alcance associado a cada ângulo (0~90)º ou (0~1.57079)rad
	double s_max;					// Alcance máximo de um lançamento
	//double y_s_max;				// Ângulo que possui alcance máximo

	double y_min, y_max, y_current;		// Ângulos mínimo, máximo e atual para intervalo de cálculo
	double error;						// Erros associados à esses intervalos
	pthread_t thread[NUM_THREADS];		// Estrutura thread para seu lançamento
	struct parameter par[NUM_THREADS];	// Parâmetros da cada thread

	// Entrada dos dados
	scanf("%c", &mode);
	if(mode == 'G'){
		scanf("%lf %lf %lf %lf %lf %lf %lf %lf", &m0, &mp, &tp, &Fp, &a, &Cd, &r, &Sf);
	}
	else if(mode == 'R'){
		scanf("%lf %lf %lf %lf %lf %lf %lf", &m0, &mp, &tp, &Fp, &a, &Cd, &Sf);
		srand(time(NULL));
		r = (double)(rand()%10+1);		// Gerando valor de raio para modo 'R'
	}
	else{
		printf("Caracter inicial inválido\n");
		return -1;
	}

	do{

		// Disparo das threads
		for(i=0; i<NUM_THREADS; i++){
			par[i].s_angle = s_angle;
			par[i].thread_id = i;
			pthread_create(&thread[i], NULL, calculate, (void*)(&par[i]));
		}

		// Aguradar as threads terminarem
		for(i=0; i<NUM_THREADS; i++){
			pthread_join(thread[i], NULL);
		}

		//Busca pelo alcance máximo
		for(i=0, s_max=0; i<NUM_ANGLES; i++){
			if(s_max < s_angle[i]){
				s_max = s_angle[i];
			}
		}

		if(Sf < s_max){	// Se for possível encontrar o alcance máximo, então calcular
			found=1;
			printf("%c", mode);
		    // Procurando o primeiro ângulo (mais baixo)
			for(i=0; !(s_angle[i] <= Sf && Sf <= s_angle[i+1]); i++);
		    y_max = (i+1)*0.01;
		    y_min = i*0.01;
			error = 1.0;
			while(fabs(error) > e){
				y_current = (y_min + y_max)/2;
				error = runge_kutta(y_current) - Sf;
				if(error > 0){
					y_max = y_current;
				}
				else{	// error <= 0
					y_min = y_current;
				}
			}
			printf(" %.5lf %.5lf", y_current, r);

			// Procurando o segundo ângulo (mais alto)
			for(i=NUM_ANGLES-1; !(s_angle[i] <= Sf && Sf <= s_angle[i-1]); i--);
		    y_max = i*0.01;
		    y_min = (i-1)*0.01;
			error = 1.0;
			while(fabs(error) > e){
				y_current = (y_min + y_max)/2;
				error = runge_kutta(y_current) - Sf;
				if(error > 0){
					y_min = y_current;
				}
				else{	// error <= 0
					y_max = y_current;
				}
			}
			printf(" %.5lf %.5lf", y_current, r);
		}
		else if(mode == 'R'){
			r *= 0.95;	//Redefine novo valor de raio
		}
	}while(mode == 'R' && !found);
	printf("\n");
	pthread_exit(NULL);
}

double gravity(double H){
	return g0*R*R/(R+H)/(R+H);
}

double density(double H){
	return p0*exp(-0.000101741*H);
}

double mass(double t){
	return ( (t <= tp) ? m0-mp*t/tp : m0-mp);
}

double force(double t){
	return ( (t <= tp) ? Fp : 0);
}

double drag(double A, double p, double V){
	return 0.5*Cd*A*p*V*V;
}




double ODE_f1(double t, double V, double y, double H){
	double D = drag(pi*r*r, density(H), V);
	return (force(t)*cos(a) - D)/mass(t) - gravity(H)*sin(y);
}

double ODE_f2(double t, double V, double y, double S, double H){
	return force(t)*sin(a)/mass(t)/V - gravity(H)*cos(y)/V;
}

double ODE_f3(double V, double y){
	return V*cos(y);
}

double ODE_f4(double V, double y){
	return V*sin(y);
}

void* calculate(void* input){
	int i;
	struct parameter *p = (struct parameter*)input;
	double *s_angle = p->s_angle;
	int id = p->thread_id;

    for(i=id; i<NUM_ANGLES; i+=NUM_THREADS){
		s_angle[i] = runge_kutta(i*0.01);
	}

	pthread_exit(NULL);
}

double runge_kutta(double y_in){
	int i;
	
	// Os 4 coeficientes do método Runge Kutta de cada EDO
	double k1_f1, k2_f1, k3_f1, k4_f1, k_f1;
	double k1_f2, k2_f2, k3_f2, k4_f2, k_f2;
	double k1_f3, k2_f3, k3_f3, k4_f3, k_f3;
	double k1_f4, k2_f4, k3_f4, k4_f4, k_f4;

	// As variáveis do problema
	double t[2], V[2], y[2], S[2], H[2];

	t[0] = 0.0;		// Tempo inicial
	V[0] = 0.1;		// Velocidade inicial
	y[0] = y_in;	// Ângulo inicial
	S[0] = 0.0;		// Deslocamento horizontal inicial
	H[0] = 0.0;		// Deslocamento vertical inicial

	for(i=0; H[i%2] >= 0; i++){

		k1_f1 = h*ODE_f1(t[i%2], V[i%2], y[i%2], H[i%2]);
		k1_f2 = h*ODE_f2(t[i%2], V[i%2], y[i%2], S[i%2], H[i%2]);
		k1_f3 = h*ODE_f3(V[i%2], y[i%2]);
		k1_f4 = h*ODE_f4(V[i%2], y[i%2]);

		k2_f1 = h*ODE_f1(t[i%2] + h/2, V[i%2] + k1_f1/2, y[i%2] + k1_f2/2, H[i%2] + k1_f4/2);
		k2_f2 = h*ODE_f2(t[i%2] + h/2, V[i%2] + k1_f1/2, y[i%2] + k1_f2/2, S[i%2] + k1_f3/2, H[i%2] + k1_f4/2);
		k2_f3 = h*ODE_f3(V[i%2] + k1_f1/2, y[i%2] + k1_f2/2);
		k2_f4 = h*ODE_f4(V[i%2] + k1_f1/2, y[i%2] + k1_f2/2);


		k3_f1 = h*ODE_f1(t[i%2] + h/2, V[i%2] + k2_f1/2, y[i%2] + k2_f2/2, H[i%2] + k2_f4/2);
		k3_f2 = h*ODE_f2(t[i%2] + h/2, V[i%2] + k2_f1/2, y[i%2] + k2_f2/2, S[i%2] + k2_f3/2, H[i%2] + k2_f4/2);
		k3_f3 = h*ODE_f3(V[i%2] + k2_f1/2, y[i%2] + k2_f2/2);
		k3_f4 = h*ODE_f4(V[i%2] + k2_f1/2, y[i%2] + k2_f2/2);

		k4_f1 = h*ODE_f1(t[i%2] + h, V[i%2] + k3_f1, y[i%2] + k3_f2, H[i%2] + k3_f4);
		k4_f2 = h*ODE_f2(t[i%2] + h, V[i%2] + k3_f1, y[i%2] + k3_f2, S[i%2] + k3_f3, H[i%2] + k3_f4);
		k4_f3 = h*ODE_f3(V[i%2] + k3_f1, y[i%2] + k3_f2);
		k4_f4 = h*ODE_f4(V[i%2] + k3_f1, y[i%2] + k3_f2);


		k_f1 = (k1_f1 + 2*k2_f1 + 2*k3_f1 + k4_f1)/6;
		k_f2 = (k1_f2 + 2*k2_f2 + 2*k3_f2 + k4_f2)/6;
		k_f3 = (k1_f3 + 2*k2_f3 + 2*k3_f3 + k4_f3)/6;
		k_f4 = (k1_f4 + 2*k2_f4 + 2*k3_f4 + k4_f4)/6;

		t[(i+1)%2] = t[i%2] + h;
		V[(i+1)%2] = V[i%2] + k_f1;
		y[(i+1)%2] = y[i%2] + k_f2;
		S[(i+1)%2] = S[i%2] + k_f3;
		H[(i+1)%2] = H[i%2] + k_f4;

		//printf("Iteraçao %d:\n", i);
		//printf("T\t\tV\t\tâ\t\tV\t\tH\n");
		//printf("%lf\t%lf\t%lf\t%lf\t%lf\n", t[(i+1)%2], V[(i+1)%2], y[(i+1)%2], S[(i+1)%2], H[(i+1)%2]);
	}

	return S[i%2];
}
